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SUMMARY 

An  analysis  of  film  condensation  on  a  vertical  fluted  tube  has  been  made 
considering  gravitational  and  surface  tension  effects  over  the  entire  fluted 
surface,  and  using  surface-oriented  coordinates.   For  the  first  time  surface 
tension  effects  are  determined,  as  they  should,  from  the  shape  of  the 
condensate-vapor  interface  rather  than  the  shape  of  the  flute. 
Two-dimensional  conduction  within  the  condensate  film  as  well  as  in  the  fluted 
tube  is  considered.   A  finite-difference  solution  of  the  highly  non-linear 
partial  differential  equation  for  the  film  thickness  is  coupled  with  a 
finite-element  solution  of  the  conduction  problem.  The  procedure  has  been 
tested  on  a  sinusoidal  flute  with  amplitude  to  pitch  ratio  O.2.     A  linear 
extrapolation,  on  a  log-log  basis,  of  our  results  shows  good  comparison  with 
experimental  data. 
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1 .   INTRODUCTION 

The  U.S.  Navy  has  a  continued  interest  in  the  reduction  of  the  size  and 
weight  of  propulsion  components  aboard  both  surface  vessels  and  submarines. 
Some  studies  are  currently  underway  in  the  Navy  to  determine  the  savings  in 
both  weight  and  volume  which  can  occur  if  the  condenser  is  designed  to  use 
corrugated  or  indented  tubing  to  enhance  heat  transfer  on  both  the  shell  and 
tube  sides.  While  a  comprehensive  research  program  has  been  underway  at  the 
Naval  Postgraduate  School  to  study  various  enhancement  techniques  for 
horizontal  condenser  tubes  [l,2],  an  attractive  alternative  aboard  submarines 
is  a  vertical  condenser  oriented  aft  of  the  steam  turbine  as  shown  in  Figure 
l(a).  With  this  arrangement,  there  would  be  a  shorter  machinery  stack  length 
compared  to  a  horizontal  condenser  mounted  underneath  the  turbine.   This  would 
allow  a  smaller  diameter  submarine  to  be  built  at  a  significant  reduction  in 
cost.  The  ultimate  aim  of  this  project  is  therefore  to  assess  the  reduction 
in  condenser  weight  and  size  that  may  be  feasible  if  the  vertical  tubes  in 
this  condenser  are  fluted  on  the  outside  (i.e.,  the  steam  side),  Figure  l(b). 

Many  methods  [3]  for  enhancing  condensation  heat  transfer  have  been 
proposed.   Among  them,  Gregorig  [k]    first  recognized  the  importance  of  surface 
tension  in  film  condensation  on  vertical  fluted  tubes.   Thereafter,  many 
experimental  studies  [5-ll]  on  vertical  fluted  surfaces  were  made  to  confirm 
his  findings.   Lustenader,  et  al.,  [5]  employed  a  fluted  vertical  tube  and 
obtained  about  four  times  larger  heat  transfer  coefficients  than  those  on  a 
vertical  smooth  tube.   Carnavos  [6]  carried  out  experiments  on  a  doubly  fluted 
vertical  tube  and  found  the  augmentation  ratio  of  condensation  heat  transfer 


rate  to  be  about  six.  Thomas  [7,8]  tested  a  vertical  tube  with  longitudinal 
rectangular  fins  and  vires.  The  augmentation  ratio  was  about  eight.   News on 
and  Hodgson  [9]  manufactured  various  condenser  tubes  of  highly  enhanced  heat 
transfer  performance.   Combs  and  his  co-workers  [10,  ll]  at  Oak  Ridge  National 
Laboratory  found  the  augmentation  ratio  to  be  around  five  for  condensation  of 
ammonia  and  refrigerants  on  vertical  fluted  tubes. 

Gregorig's  theoretical  model  has  also  been  improved  upon  but  much  still 
remains  to  be  done.   Edwards,  et  al.,  [12]  proposed  a  condensation  model  on  a 
heat  transfer  surface  of  triangular  fins  on  the  assumption  of  liquid  film 
attaching  to  the  tip  of  fin  with  a  contact  angle.  The  effect  of  a  locally 
thin  condensate  film  on  the  side  of  the  fin  (Figure  2)  was  not  considered  by 
them  as  well  as  by  Fujii  and  Honda  [13] .   The  latter,  however,  did  consider 
two-dimensional  conduction  within  the  condensate  film  in  the  trough  region 
(Figure  2)  and  within  the  fin  wall.  Mori,  et  al.,  [lU]  considered  the  effect 
of  a  thin  condensate  film  on  the  side  of  the  fin  but  neglected  the  variation 
of  its  thickness  in  the  vertical  direction.   They  also  considered  only 
one-dimensional  conduction  within  the  fin  and  neglected  any  conduction  within 
the  condensate  film  in  the  trough  region.   Hirasawa,  et  al.,  [15]  improved 
upon  [lU]  in  that  they  included  the  variation  of  a  thin  condensate  film  in  the 
vertical  direction  but  neglected  conduction  within  the  fin  and  the  film 
completely.   Panchal  and  Bell  [l6,  17]  also  neglected  conduction  within  the 
fin  and  the  film  while  analysing  a  sinusoidal  fluted  tube,  but  later  found 
that  two-dimensional  conduction  is  important  within  the  fin  and  the  film  for  a 
triangular  fin  [l8] .  A  recent  empirical  analysis  by  Barnes  and  Rohsenow  [19] , 
based  largely  on  [l6,  20],  reports  an  augmentation  ratio  of  about  fifteen  for 
condensation  of  steam  on  a  fluted  surface  while  most  earlier  studies 


considered  condensation  of  a  refrigerant  and  found  much  smaller  augmentation 
ratios.  With  present  day  knowledge,  it  is  difficult  to  ascertain  whether  this 
discrepancy  is  due  to  different  fluids  or  due  to  a  questionable  analysis. 

All  theoretical  analyses  discussed  above  break  up  the  fluted  surface  into 
basically  two  parts.   In  the  portion  near  the  crest,  gravity  is  neglected  in 
comparison  to  the  surface  tension  effect,  while  in  the  portion  near  the 
trough,  the  reverse  is  done.   The  patching  between  the  two  regions  is  carried 
out  at  a  point  that  is  selected  quite  arbitrarily  at  times  [ 15-18] .   This 
isolation  of  the  two  important  effects,  gravity  and  surface  tension,  is 
justified  on  the  basis  that  the  condensate  film  is  thick  in  the  trough  region 
and  thin  over  the  crest.   This  is,  however,  not  true  over  the  initial  portion 
of  the  tube  length  in  the  vertical  direction  nor  during  the  initial  portion  of 
the  tube  just  below  a  condensate  drainage  skirt  (Figure  1(b)).   These  initial 
portions  may  be  a  few  pitches  in  length.  A  recent  analysis  by  Stack  and 
Merkle  [21 ]  does  attempt  to  solve  the  complete  equation  with  both 
gravitational  and  surface  tension  effects  included  over  the  entire  flute  but 
has  three  major  drawbacks.   First,  when  the  condensate  film  in  the  trough 
region  has  thickened,  the  analysis  treats  it  in  the  boundary-layer-sense,  that 
is,  it  neglects  all  velocity  gradients  except  those  in  the  direction  normal  to 
the  fluted  surface.  Second,  the  analysis  is  restricted  to  impracticably  low 
values  of  amplitude  to  pitch  ratio  of  the  flute  (0.02  and  0.0U)  owing  to  the 
use  of  a  Cartesian  coordinate  system  rather  than  a  surface-oriented  coordinate 
system.   Third,  the  analysis  does  not  consider  any  heat  transfer  effects. 
Their  analysis  is  only  confined  to  finding  the  condensate-vapor  interface. 
Fujii  and  Honda  [13]  also  considered  the  entire  thin  film  as  one  piece  over  an 
initial  length  of  the  tube  (about  1/8  of  the  pitch)  but  neglected  the 
important  surface  tension  effect. 


A  major  deficiency  of  all  the  above  theoretical  analyses  except  [21 ]  lies 
in  the  way  surface  tension  effect  is  determined.   Contrary  to  the  general 
belief,  this  effect  does  not  depend  upon  the  curvature  of  the  condensate-vapor 
interface.   Instead  it  depends  upon  the  variation  of  this  curvature  along  the 
interface.   However,  since  the  location  of  this  interface  is  unknown  a  priori, 
many  analyses  [12,  13]  simply  determine  it  on  the  basis  of  the  known  flute 
shape  owing  to  the  argument  that  in  the  crest  region  where  the  surface  tension 
effect  is  important,  the  film  is  very  thin.   Such  analyses  cannot  be  applied 
at  all  to  triangular  or  rectangular  fins  since  their  curvature  as  well  as  the 
variation  of  curvature  along  the  flute  is  zero.  The  same  is,  however,  not 
true  for  the  condensate-vapor  interface.  Even  for  sinusoidal  flutes,  there 
are  large  differences  between  the  curvature  and  its  variation  along  the  curved 
surface  for  the  given  flute  and  the  actual  interface  (see  Fig.  10  and  its 
discussion  in  Section  5)»   Moreover  those  analyses  that  do  claim  to  determine 
the  surface  tension  effect  based  on  the  actual  shape  of  the  condensate-vapor 
interface  are  also  questionable  in  their  claim  as  shown  in  Appendix  A. 

Another  problem  with  all  earlier  theoretical  analyses  except  [13]  is  that 
conduction  within  the  condensate  film  and  the  tube  wall  is  either  completely 
neglected  or  considered  to  be  at  most  one-dimensional.   Panchal  and  Bell  [l8] 
point  out  clearly  that  two-dimensional  conduction  should  be  considered  within 
both  the  condensate  film  and  the  tube  wall. 

In  an  attempt  to  take  care  of  these  deficiencies,  an  analysis  has  been 
developed  during  this  study  that  solves  the  complete  equation  with 
gravitational  as  well  as  surface  tension  effects  included  over  the  entire 
fluted  surface  using  the  surface-oriented  coordinate  system.   In  addition,  for 
the  first  time,  surface  tension  effects  are  determined  as  they  should,  from 


the  actual  shape  of  the  condensate-vapor  interface  rather  than  the  shape  of 
the  flute.   Two-dimensional  conduction  within  the  condensate  film  as  well  as 
the  fluted  tube  is  also  included.   A  finite-difference  solution  of  the  highly 
non-linear  partial  differential  equation  for  the  condensate  film  thickness  is 
coupled  with  a  finite-element  solution  of  the  conduction  problem.   Details  of 
the  analysis  follow  in  subsequent  sections. 


2.   ANALYSIS 

Considering  the  condensate  as  a  viscous,  incompressible  Newtonian  fluid, 
and  setting  up  a  curvilinear  orthogonal  coordinate  system  (x.,  x  ,  z)  as  shown 
in  Fig.  3,  the  continuity  equation  yields  [22] 
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where  (u  ,  u  ,  w)  are  the  velocity  components  in  the  (x  ,  x  ,  z)  directions 
respectively,  and  R(x  )  is  the  radius  of  curvature  of  the  fluted  tube.  The 
momentum  equations  are  (cf.  Ref.  22,  p.  68) 


X]_-direction 


3u, 


u. 


R+x   1  3x 


3u, 


u 


2  3x, 


w 


3u1 

+ 

U1U2 

+ 

1     R        3p 

3z 

R+x2 

p  R+x2    3x1 

=   V 


R2         3S 

+ 

3x2 

+ 

3z2 

+ 

1 

R+x2 

3ux 
3x2 

Ul 

(R+x2)2    3x2 

(R+x2)2 

2R 


3u, 


(R+x2)2  8xl 


R 


dR 


u. 


(R+x2)3dxl  2 


Rx, 


dR 


3u. 

"3Z 


(R+x2)3  dxl 


(2.2) 


X2~direction: 


R 


3u, 


3u, 


u. 


R  +   x        1   3x. 


u 


2    3x, 


w 


c 

~3z~ 


u. 


R  +  x. 


+       — 


1  _3p_ 

P    3x, 


=      v 


R2              *\ 

+ 

s2 
3  u2 

3x2 

+ 

32u2 
3z2 

+ 

1          ^2 

(R   +  x2)2    3x2 

R  +   x2    3x2 

u. 


2R 


3u, 


(r+x2v 


(R+x2)2   8X1 


R  dR  +  RX2        dR      ^2 

(R+x2)3  dXl  Ul  (R+x2)3   dxl    3X1 


(2.3) 


z-(vertical)    direction 
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where  p  and  v  are  the  density  and  kinematic  viscosity  of  the  condensate,  p  is 
the  pressure  in  the  condensate  film,  g  is  the  acceleration  due  to  gravity,  and 
pv  is  the  density  of  the  vapor. 

In  order  to  simplify  these  equations,  we  make  the  usual  assumptions  that 
inertia  terms  are  negligible  compared  to  other  terms,  and  that  3/ 3x  »  3/ 3x 


or  3/3z,  and  vu  <<  u.  or  w.   This  simplifies  the  momentum  equations  (2.2)  - 
(2.U)  to 
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where  u  (=pv)  is  the  dynamic  viscosity  of  the  condensate.   Integration  of 
(2.7)  with  boundary  conditions 
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Here  S  is  the  (known)  shear  stress  (in  z-direction)  on  the  condensate-vapor 
interface,  and  5(x  ,z)  is  the  condensate  film  thickness.   Similarly, 
integration  of  (2.5)  with  boundary  conditions 
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Here  S  is  the  (known)  shear  stress  (in  x  -direction)  on  the  condensate-vapor 
interface. 

It  is  advantageous  to  use  the  integral  form  of  the  continuity  equation 
rather  than  the  differential  form  (eqn.  (2.1)).   We  therefore  consider  the 
control  volume  of  Fig.  H,  and  write  the  continuity  equation  as 
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where  kf  is  the  thermal  conductivity  of  the  fluid,  and  hf„  is  the  latent  heat 
of  vaporization.  The  right  side  of  eqn.  (2.12)  represents  the  rate  at  which 
condensation  of  vapor  takes  place.  The  temperature  gradient  3T/ 3x  can  be 


11 


approximated  by  (Ts  -  Tw)/5,  where  Ts  is  the  saturation  temperature  at  which 

condensation  takes  place,  and  Tw  (x  ,z)  is  the  wall  temperature  at  x  =  0. 

This  approximation  basically  assumes  one-dimensional  conduction  in  the 

condensate  film  but  since  3/3x  >>  3/3x.  or  3/3z,  and  we  are  interested  only 

in  (3T/3x  )  at  present,  it  is  justifiable  to  make  this  approximation. 
2  6 

Substituting  for  w  and  u  from  eqns.  (2.9)  and  (2.11)  into  eqn.  (2.12), 
and  integrating,  we  get 
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The  pressure  p  can  be  related  to  the  surface  tension  a  and  radius  of  curvature 
R-j_  of  the  condensate  vapor  interface  by 


p  =  ps  ±  a/Ri   , 


(2.1U) 


where  ps  is  the  saturation  pressure  of  the  vapor.  The  positive  sign  holds  for 
0  <  x  <  xp/2>  and  the  negative  sign  for  Xp/2  <  x   <  Xp,  where  Xp  is  the 
length  (curve  DE  in  Fig.  3)  along  the  fluted  surface  (in  x  -direction)  from 
the  crest  to  the  trough  (over  half  the  pitch).  We  are  basically  interested  in 
analyzing  flute  shapes  that  are  symmetric  about  the  crest  and  trough,  and 
therefore  take  advantage  of  the  symmetry  by  considering  only  half  of  the  flute 
pitch. 
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Assuming  a   to  be  constant,  the  pressure  gradient  may  be  calculated  from 
eqn.  (2.lU) 

air    -    *  °  &(k)'  (2-15) 

1  1  v  1  ' 

It  is  thus  the  variation  of  curvature  of  the  condensate-vapor  interface  that 
is  of  significance,  not  the  curvature  itself.   In  general,  this  can  be  quite 
different  from  d(l/R)/dx  ,  a  fact  neglected  by  many  previous  analyses. 
Relations  for  finding  d(l/Rj_)/dx  are  given  in  Appendix  A. 

Equation  (2.13)  is  a  partial  differential  equation  for  the  condensate 
film  thickness  <5(x  ,z).   It  requires  the  prior  knowledge  of  shear  stresses  S, 
and  S  that  come  from  vapor  dynamics  and  of  Tw  (x  ,z)  that  comes  from  heat 
conduction  analysis  for  the  fluted  tube  and  condensate  film.   Eqn.  (2.13) 
involves  only  first  order  derivatives  w.r.t.  z  but  fourth-order  derivatives 
w.r.t.  x.  owing  to  the  presence  of  dp/dx   (cf.  eqn.  (2.15)  and  (A-15)  or 
(A-l8)).   For  flute  shapes  that  are  symmetric  about  the  crest  and  trough,  the 
boundary  conditions  for  the  solution  of  eqn.  (2.13)  are 


85       335 


-   =   0  at  x  =  0,  and  at  x  =  x  for  all  z;       (2.l6a) 
3  1  1    p 


9X1      ax3 


5(Xl,0)  =  0. 

While  the  boundary  condition  at  z  =  0  is  correct,  it  is  not  practical  to  start 
the  integration  of  eqn.  (2.13)  from  z  =  0.   Therefore,  following  Stack  and 
Merkle  [2l]  ,  we  replace  the  boundary  condition  at  z  =  0  by 

5(Xl,z0)  =  6Q    ,  (2.16b) 
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where  6     is  the  film  thickness  at  the  initial  station  z   .  and  is  found  by  the 
0  0  - 

classical  two-dimensional  Nusselt  solution 

1/4 


4yk.z  (T  -T  J 

*   °l    °     "°°j  (2.17) 

p(p-p„)  gnfg    ( 


v 
Here  Tw00   =  Tw  (0,z0)   »  Tw  (0,0). 

For  non-dimensionalization,  we  take  the  length  L  of  the  vertical  fluted 

tube  in  the  z-direction  and  length  Xp  and   6y  as  the  characteristic  lengths  in 

the  x   ,  and  x  -directions  respectively.     Here   5y  is  related  to  L  by  the 

classical  Nusselt  relation   (2.17)   in  exactly  the  same  manner  as   6     is  related 

to  z„ .     Thus ,  we  let 
0 

A  =  6/6      ,        X  =  x,/x      ,        Z  =  z/L     .  (2.18) 

v  *■     P 

With  this  non-dimensionalization,  eqn.    (2.13)   and  boundary  conditions   (2.16)     . 
can  be  written  as 


3A3        _     3Q  _       8A2  _        8A2  A2/3E1  3E3  \  3°2   ,      ,01Q. 

+  D.  —    +    E^     —-    +      E_    ~=-    +     a  hv-    +      w      =    mr        (2-19) 


3Z  1  3X  1     3X  3      3Z  \  3X  3Z     J  4A 


3 
|4    =    -^— I      =0       at  X  =  0,1     for  all  Z,  (2.20a) 

3X  3XJ 


1/4 
A(X,ZQ)    =  ZQ  ,  (2.20b) 


where 
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zo  =  VL     ' 


Dl  = 


3a  L 

2      ' 

g(p~P    )    *     5 
3   K    *v       p     y 


D2  = 


T         T 
s  -    w 

T     -  T  ' 

S  WOO 


3IS-, 


Ei  = 


2g  (p-p   )x  6 


3S. 


E„   = 


2g(P-Pv)6y      ' 


Q  =    Re 


Pc  A/2  +  3  D    A2/4  -   (PcA  +  3§/2D  )   In    (1+D  A/I^) 


F      , 


D3   =     6y/xp     , 


Fq  =    R/xp    / 


F  =     ± 


and 


d_ 
dX 


(k)'^AH- 


-  f or  0   <  X  <  1/2 
+  for  1/2  <  X  <  1 


Pn  =     Ri/xp 

Equation   (2.19)    is  a  highly  non-linear  partial  differential  equation.      It 
is  solved  numerically  by  the  finite-difference  technique  as  detailed  below. 
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3.   FINITE-DIFFERENCE  METHOD 

Equation  (2.19)  is  parabolic  in  Z  so  that  a  forward  marching  scheme  in 
the  Z-direction  can  be  used.   Thus  the  choice  of  6  will  affect  only  the 
region  near  z  ,  and  its  effect  will  die  out  as  z  increases.   This  was  indeed 
found  to  be  true.  Equation  (2.19)  as  written  is  in  conservative  form.   An 
equivalent  but  non-conservative  form  can  be  obtained  by  expanding  the 
derivatives  in  (2.19)  and  dividing  by  3 A2.   In  that  non-conservative  form,  the 
equation  is  still  parabolic  in  Z  but  the  mass  is  not  identically  conserved 
when  the  equation  is  integrated  numerically  [2l] . 

We  divide  the  interval  0  <  X  <  1  into  n  equal  parts.  Taking  5X  and  <5Z  as 
the  step  sizes  in  X  and  Z-directions  respectively,  and  using  backward 
differencing  in  Z  and  mixed  differencing  in  X-direction,  we  write  the 
finite-difference  form  of  eqn.  (2.19)  as 


A3i,k-  A3i,k-1    n   gQi+l,k  +  (l-2g)Qi,k-(l-g)Qi-l,k 

Vsz)  +  Di  Tsx) 


3A2i+lsk+(l-23)A2ijk-(l-3)A2i_l5k        A2^  _  A2^ 

Ei  Tsx]  +  s3      Vsz) 


2     /  8E1     9E3  \        3D2 

1  ,K 
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where  the  subscripts  i  and  k  represent  the  location  in  X  and  Z  directions 
respectively  (cf.  Fig.  5)»  and  3  is  a  number  that  is  selected  between  0  and  1; 
3=0  corresponds  to  backward  differencing,  3  =  1  to  forward  differencing,  and 
3  =  1/2  to  central  differencing  in  X.   Since  eqn.  (2.20b)  gives  A's  for  all  X 
(i.e.,  for  all  i)  at  Z   (say  k  =  l) ,  we  can  march  forward  in  the  Z-direction. 
In  other  words,  A's  at  (k-l)  for  all  i  are  known  in  eqn.  (3.1).   Thus  (3.1) 
leads  to  a  tridiagonal  set  of  non-linear  algebraic  equations  to  be  solved  for 
(n+l)  values  of  A's  at  each  location  k  of  forward  march  in  the  Z-direction. 
This  non-linear  set  is  solved  by  linearization  and  successive  iteration.   Let 
us  write  eqn.  (3.1)  in  the  form 

A.(A.  _  Jm  +  B.(A.  \m     +  C.(A.+1  ,\m  =  R.       ,    (3.2) 
i\  l-l ,k/       i\  i,k/       i\  l+l ,k;       l 

where  (  )m  represents  the  values  at  the  (current)  mth  iteration,  and  the 
coefficients  A^ ,  B-j_ ,  C-j_ ,  and  Rj_  depend  upon  the  method  of  linearization. 

For  linearization  three  methods  were  tried.   In  the  first  method,  we 
simply  let 


.    m       /   4  i   m-1         m 

(»i.*)    ■  (4i;i)     K,J    ■  j- 1.2.3.       o.3) 


where  (  )m_1  represents  the  known  values  at  the  (previous)  (m-l)th  iteration. 
The  iteration  is  started  by  taking  A-j_  jj  =  A^  ^_]_  for  all  i.  With  this 
linearization  the  coefficients  in  eqn.  (3.2)  are 


IT 


*!-'-!§  (1-B) 


iK-i>kr  +  Ei(vi,k) 


ra-1 


,    i  =  l,2,...,(n-l)        (3.Ua) 


Bi  =  (Ai,k) 


m-l 


K.J"         +E3+«ZV 


3E,         3E, 

+ 


i,kJ 


♦     §     d-23) 


'iKJ 


m":         +      E,    (  A.    v  ^ 
1  \    i,k 


,   i  =  0,1.. .  ,n     (3.Ub) 


Ci    "      5X   3 


Dl(G 


i+l,k 


m"x  +  MW1 


i   =  l,2,...,(n-l)    (3.Uc) 


o  p  3D   (6Z) 

Ri  =       Ai,k-1  +  E3   Ai,k-1  +     7T        ,m-i        •   *  =  0,1.-., n  (3.Ud) 

MAi,k) 


A     =     -^ 
n  SX 


D.(G     1    .  r*1  +  E_    (1-2  6)    (A     .    . 
1\    n-l,k/  1  \   n-1,] 


m-l 


(3.Ue) 


r     =     i£ 

0  5X 


i.(°i.k)      -  V1-2"  K.J 


\m-l 


(3.Uf) 


where 


G  =  R 


R   /2  +  3D_A/U 
c  3 


-( 


Rc  +   2  Dl    A 


In  (     1   +  D   A/R  ) 


F  =   Q/A        ,  (3.Ug) 
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and  the  boundary  conditions  in  (2.20a)  have  been  taken  into  account.   This 
method  works  but  has  two  drawbacks.   Its  convergence  is  slow,  and  it  requires 
heavy  under-relaxation  for  convergence.   The  latter  implies  that  it  is 
necessary  to  under-relax  the  values  of  A  after  every  iteration  according  to 


m-1 


(A)ra  =  (A)m-X  +  X 


A)IU-  (A) 


m-1 


(3.5)* 


where  the  relaxation  factor  X  <  1.   In  fact,  it  was  found  that  X  decreases 
from  about  0.5  at  Z  to  very  small  values  as  Z  increases.   Typically  X  <  0.1 
for  Z  >  0.02. 

In  the  second  method,  we  used  Taylor  series  expansion  in  order  to 
linearize  terms  containing  A2  and  A3.  Thus  we  let 


-i    hi       in*  m— 1       m       ,  ~ 

(4, J    -  3KJ    K,J  -  *Kk) 


m-1 


(3.6a) 


and 


m-1 


(3.6b) 


With  this  linearization  the  coefficients  in  eqn.  (3.2)  are 


A.  =  -  ||  (1-6) 

1         OA 


Dl(Gi-l,k) 


m-1 


+  2  E, 


iK-i,k) 


m-1 


,  i=l,2,...,(n-l)   (3.7a) 


*  This  equation  is  to  be  treated  in  the  context  of  an  equivalent  FORTRAN 
statement . 
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=  (vr-1  3(4i)kr1  +  2E3  +  2(5z)(-3r  +  ir) 


l,k-J 


♦  f  (1-2B) 


,111-1 


i(  Gi>k)       +  2  Ei(  \*) 


m-1 


,     i  =  0,1,. . .  ,n  (3.7"b) 


r    -    ^    8 


Dl(Gi+l,k) 


m-1 


2  Ei     (   Ai 


1    \     i+l,k 


m-1 


i   =   1,2,. ..(n-1)      (3.7c) 


3D    (6Z)  -       m-1 

R.    =  A.    .     .    +  E0   AT  .     -    + r  +  2   A.    . 

l  i,k-l  3     i,k-l       ,,/A        \m-l  V   i,k' 


*K.k) 


+   El   6X 


m-1 


m-1 


"(W)     +(1-26)(^,J     -  (1-6)  (4-i,J 


m-1 


3E  3E 

E3  +  6Z    l^T^ 


/  3E  3E_  \        "I       p  m-1 

(^T+^).    ,      (Ai,k)         •   i=l>2,...,(n-l)        (3.7d) 


i,kJ 


a   =   -i* 

n  6X 


D.(G     .    v)m_1   +  2E.    (1-23)   (A     .    .    ) 
1\    n-1  ,k/  1  \    n-l,k  / 


m-1 


(3.7e) 


L0  5X 


'it5!,,)1""1      -     ^l-^2  *»(*!,, 


m-1 


(3.7f) 
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Ro  '  4o,k-i  +  E3  'o.k-1 


3DP(,SZ)  A7  ,  9    1 


m-l 


MvX 


«,r 


2  A 


0,k 


)m-1  +  B^  ||  (1-23)  +  E3  +  5Z 


3E. 

~3X~ 


3E, 

Iz" 


0,k-i 


(3.7g) 


3D2(5Z) 


5Z 


3  2 

n    n,k-l    3  n,k-l   ,,  /  A    \m-l   "1  6X 


-  -  E,  -7^  (1-2  3)  (a  ,  , 
\  n-1 ,1 


m-l 


\    n,k/ 


i-i  r 


(tl  Mv*)""1 


+  E1  ||  (1-23)  +*E3  +  5Z 


3E 


3E. 

~3z~ 


n,k  _ 


(3.7h) 


where  G  is  again  given  by  (3.^+g)  and  the  boundary  conditions  in  (2.20a)  have 
been  accounted  for. 

It  is  of  interest  to  note  from  eqns .  (3.7)  that  the  Taylor  series 
expansion  was  not  used  to  linearize  the  right  side  of  eqn.  (3.1)  as  well  as 
the  terms  involving  Q  on  the  left  side  of  (3.1).  While  such  a  linearization 
of  the  right  side  of  (3.1)  can  be  carried  out  without  much  difficulty,  it  was 
not  found  to  be  beneficial.  The  unnecessary  complication  was  therefore 
avoided.   The  Taylor  series  linearization  of  terms  involving  Q  in  (3.1)  was 
also  tried  but  was  found  not  to  work  at  all.   This  method  was  found  to  be  the 
best  of  the  three.   It  allows  a  larger  step  size  ( 6Z)  in  the  Z-direction,  and 
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a  larger  relaxation  factor  X  than  the  first  method.   Overall,  it  is  better 

than  the  first  method  by  a  factor  of  at  least  twenty,  and  was  therefore 

adopted.   The  third  method  consisted  in  the  use  of  Wewton-Raphson  method  for 

solution  of  the  non-linear  equations.  A  general  outline  of  this  method  is 
given  in  Appendix  B.   Though  this  method  does  not  require  any  relaxation 

factor,  it  produced  absurd  results  (ATs  oscilling  with  X  after  a  few  step 

sizes  in  the  Z-direction) ,  and  was  therefore  abandoned. 
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k.      COMPUTATIONAL  DETAILS 

Before  the  set  of  equations  (3.2)  resulting  from  eqn.(3.l)  can  "be  solved, 
we  need  to  specify  values  of  the  dimensionless  parameters  E  ,  E  and  D  .  This 
requires  the  prior  knowledge  of  shear  stresses  S  and  S  ,  and  the  temperature 
Tw(x  ,  z)  of  the  condensate-wall  interface.   The  computer  code  does  have  the 
provision  for  specification  of  shear  stresses  S  and  S  but  for  the  present 
results  both  these  stresses  were  set  to  zero.  This  is  "because  vapor  dynamics 
that  determines  these  stresses  is  beyond  the  scope  of  the  present  work.   The 
computer  code,  however,  does  calculate  Tw(x  ,z)  by  considering  two-dimensional 
conduction  within  the  fluted  tube  as  well  as  the  condensate  film.  For  this, 
the  Laplace  equation 

k(i!|   +if|\  .    0  (U.1) 

V  3x     3y  ' 

is  solved  subject  to  the  boundary  conditions 

8T 

7-  =  0  at  x  =  0,  P/2  where  P  =  pitch  of  the  flute, 

oX 

T  =  Ts  at  the  vapor-condensate  interface,  (^.2) 

and    T  =  Tc  at  the  coolant-tube  interface , 
where  k  is  the  thermal  conductivity,  and  Tc  is  the  coolant  temperature  assumed 
to  vary  linearly  from  the  coolant  inlet  to  exit  temperature.  The  last 
boundary  condition  in  (U.2)  needs  modification  due  to  a  film  resistance  on  the 
coolant  side  but  such  a  specification  either  requires  an  ad-hoc  value  of  the 
film  resistance  or  a  consideration  of  coolant  dynamics  that  was  also 
considered  beyond  the  scope  of  the  present  analysis.   Solution  of  eqns .  (U.l) 
and  (H.2)  was  obtained  by  a  finite  element  method  in  which  the  region  0ABC0 
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(Fig. 3)  was  divided  into  several  linear  triangular  elements.  Details  for  this 
can  be  found  in  any  text  on  the  finite  element  method.  Some  care  is,  however, 
required  since  the  thermal  conductivity  of  the  fluid  in  region  ABEDA  is  vastly 
different  from  that  of  the  tube  material  in  region  DECOD  (Fig. 3).  The  finite 
element  method  was  preferred  over  the  finite-difference  method  owing  to  the 
irregular  shape  of  surface  AB  (Fig.  3).  From  this  solution,  we  get  the  values 
of  Tw(x  ,z)  on  the  surface  DE  (Fig.3). 

Since  the  solution  of  eqn.  (3.1)  requires  Tv(x  ,z)  which  comes  from  the 
solution  of  (U.l)  for  which  knowledge  of  6-values  is  essential  (for  locating 
surface  AB  in  Fig.3),  the  two  solutions  are  coupled.   Also,  since  the  solution 
of  (3.1)  is  found  iteratively,  it  implies  that  eqn.(U.l)  should  be  solved  at 
every  iteration.   This  places  a  rather  prohibitive  demand  on  the  computer 
time.   However,  since  (6Z)  is  small,  of  the  order  of  10" 5,  changes  in  5  are 
small  at  every  iteration  with  the  result  that  it  is  not  really  necessary  to 
solve  equation  (*+.l)  at  every  iteration.   Thus  (U.l)  was  solved  only  once  per 
step  in  the  Z-direction.  This  saves  considerable  computer  time  since  the 
finite  element  solution  of  (U.l)  requires  almost  as  much  computer  time  as 
nearly  one  hundred  iterations  for  solution  of  (3«l).   Some  sample  runs  were 
initially  made  to  confirm  that  the  error  in  6-values  due  to  this  time  saving 
feature  was  really  negligible. 

U.l  Condensate-Vapor  Interface  Profile 

As  pointed  out  in  Appendix  A,  one  method  for  finding  d(l/Rj_)/dx  is  by 
use  of  equation  (A. 8)  which  involves  derivatives  of  f-j_(x),  where  fj_(x) 
describes  the  condensate-vapor  interface  at  any  Z.   Several  methods  were  used 
to  find  these  derivatives  but  none  of  them  proved  worthwhile  since  fj_(x)  is 
known  at  unequally-spaced  values  of  x.  The  methods  are  listed  below  for  the 
sake  of  completeness. 
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i)   Cubic  splines  were  fitted  in  a  least  squares  manner  through  f-j_  values 
to  find  fj_  and  fj_  .   Then  cubic  splines  were  fitted  through  the  previously 
computed  f-j_  values  to  find  fj    .   A  slight  variation  of  this  method  was  to 
fit  cubic  splines  repeatedly  through  f -•_ ,  f^  and  f  ■•_  values  to  find  f ± 
Both  these  methods  gave  different  values  for  fj_  and  fj    pointing  to  their 
futility. 

ii)   Parabolas  were  fitted  through  three  consecutive  points  in  succession 
on  the  condensate-vapor  interface,  and  R^  and  d(l/R-j_)/dx  were  computed  at  the 
point  central  to  each  set  of  three  points.  This  also  led  to  inconsistent 
results . 

iii)  An  attempt  to  fit  a  truncated  Fourier  series  (having  3-^  terms) 
through  fj_  values  was  made.   This  was  attempted  since  we  were  working  with  a 
sinusoidal  flute.   However,  since  the  fit  itself  was  poor,  no  attempt  was  made 
to  calculate  the  derivatives.   A  better  fit*  (Fig.  6)  was  obtained  by  fitting 
(a  +  a  cos  2ttx/P)  through  all  (n+l)  values  of  f-j_,  where  the  coefficients  a 
and  a  were  found  for  a  least  squares  fit.   However,  even  this  fit  was 
considered  inadequate  to  find  d(l/R-j_)/dx  .  A  slight  variation  of  this  was  to 
fit  (a  +  a  cos  2ttx/P)  through  the  first  n/2  values  of  fj_  (for  x  <  P/2)  and  a 
similar  curve  through  the  last  n/2  values.   This  produced  a  very  good  fit  on 
visual  inspection  (Fig. 7).   However,  using  it  to  find  d(l/R-j_)/dx  leads  to 
wiggles  in  6 -values.   Solution  of  equations  (3.1 )  and  (U.l)  starting  from  ZQ 
leads  to  film  shape  rising  suddenly  in  the  trough  region  (Fig. 8)  after  a  few 
(5Z)  steps,  thus  violating  the  boundary  conditions  there.   Starting  from 
Z  =  0.07,  with  results  at  Z  =  0.07  found  by  better  methods,  the  solution  leads 
to  such  large  wiggles  that  even  negative  values  of  6  are  found  only  a  few 
steps  downstream.   In  an  effort  to  improve  upon  this  situation,  values  of  6 
were  computed  to  correspond  to  the  fit  in  Fig.  7  before  iterating.  This, 
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however,  slowed  down  the  convergence  so  much  that  even  two  hundred  iterations 
were  insufficient  for  the  same  tolerance.   In  fact,  it  appeared  as  if 
convergence  would  never  take  place.   Nevertheless,  this  exercise  did  point  out 
the  reason  for  the  failure  of  the  above  technique,  since  the  differences 
between  the  original  5-values  and  those  conforming  to  the  good  ( ! )  fit  in 
Fig.  T  were  found  to  be  large,  as  much  as  20$. 

iv)   Efforts  to  use  lagrangian  interpolation  to  compute  f-j_  values  at 
equi-distant  x  values,  and  then  use  techniques  described  in  Section  U.2  to 
find  f-j_  ,  etc.  also  produced  absurd  results.  Use  of  equation  (A. 8)  to 
compute  d(l/Rj_)/dx  was  therefore  abandoned. 

It  may  be  mentioned  that  while  using  eqn.  (1.8)  to  find 
d(l/R-i_)/dx.  ,  a  futile  attempt  was  also  made  to  solve  the  partial  differential 
equation  (2.13)  directly  in  terms  of  63  rather  than  6. 

U.2  Derivatives  of  Condensate  Film  Thickness 

Use  of  either  eqn.  (A. 15)  or  (A.l8)  to  find  d(l/Rj_)/dx  requires  the 
determination  of  first  three  derivatives  of  6  w.r.t.  x  .  Fortunately  5  is 
known  at  equidistant  values  of  x  .   It  is  therefore  easier  to  find  these  der- 
ivatives than  those  of  fj_  w.r.t.  x.   Several  methods  for  this  determination 
were  also  tried  with  varying  degrees  of  success. 

i)   Derivatives  of  delta  were  found  by  fitting  a  cosine  to  the  vapor- 
condensate  interface  in  a  manner  analogous  to  that  described  in  Sec.  *+.l,  but 
it  resulted  in  absolutely  useless  results. 

ii)  Another  fruitless  effort  was  to  fit  least  squares  cubic  splines, 
repeatedly  or  otherwise ,  through  6-values  so  as  to  find  6  ,  5  and  5 
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iii)  The  first  derivative  of  5  w.r.t.  x  was  found  by  solving  the  following 
set  of  linear  equations  [23,  p.56]  : 


5i-l  +  k<±   +  5i+l  =^(5i+l  "  Vl  )   -  i  =  l^,...,(n-l), 


(U.3) 


5=0    ,    5  =  0   ,   5x.  =  x  (6X), 
o  n  1    p 


and  the  second  derivative,  5",  by  solving 


If  II        M 


6.  ,  +  105.  +  6, 


l-l   __„    ^i+1 


K)! 


(  5i-l  "  26i  +  Vl)'   i=l,2,...,(n-l), 


55"  +  s"   =  t^-,o  (5,-5  I  ,  (k.h) 


1   K 


\2  \  1   o  /  ' 


m        11     jo     / 

Vx^n— )2(Vl- 
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Though  these  sets  of  linear  equations  are  tridiagonal  and  can  be  solved 

easily,  they  do  represent  an  increase  in  computer  time.   Besides  the  method 

yields  rather  poor  results. 

iv)   A  method  that  gives  better  results  than  those  obtained  by  above 

methods  is  based  on  passing  least  squares  quartics  through  seven  consecutive 

points  at  a  time.   Leaving  the  details  to  the  reader  (see  also  [2U] ,  p.^92), 

t    ii      fii 

we  present  the  equations  for  finding  6j_,  Sj_  and  6-j_   for  i  =  o,  1,...,  n. 

These  equations  make  use  of  the  boundary  conditions  (2.l6a),  and  the  fact  that 
5  is  known  at  equidistant  x,  values. 

Relations  for  5',  accurate  to  0(5x  J  ,  are 

6't  ■  (226i-3-  676i-2"  58Vi+  58W  67<W-  ^Vs)'  2526xi  • 

i  =  3,^,.  •  •  ,(n-3) , 

«'    =  0      ,      5*    =  0S 

o  n 


61   =  \'586o  "  6T61+  8062  +  6T53  "  22\)/   2^26\    » 
52   =  ("6T5o   "  3661+   5863   +  6lSk  "  2265)/   2526xi    » 


6*        =  (226     _   -  676     ,    -  586     _   +  365     .    +  676    )/   2525x,    , 
n-2        \        n-5  n-4  n-3  n-1  n/  1 


(U.5) 


and 


6*        =  (225  -  675     .   -  806     _  +  67  5     .    +  58  5    )/   252 6x. 

n-1       \        n-k  n-3  n-2  n-1  n/  1 


-28- 


Relations    for   6",    accurate   to   01  5x    r,   are 

5i  =  ("136i-3+  6T5i-2"  196i-l"  T05i  "  19W  6T5i+2  "  135i+3)7   1326X1    ' 

i  =  3,!+,...,(n-3), 

6"  =  (-265^+  13H60  -  386n    -  706    )/   1325x2    , 
o        \  3  2  1  0/  1 


61  =  (-135U+  6T63  "  3262  "  351"  196o)/   1325xi    ' 


6^  =  (-1355+  675^  -  19<53  -  T052-  326^  67  6q j /   1326x2    ,  (U.6) 


6*'   .   =  (-136     «.+  676     ,,-  196     ,-  706     _-326     .+  676  )/   1326x2    , 
n-2       \  n-5  n-h  n-3  n-2  n-1  n/  1 


5"  .    =  (-136     ,+  675     ,-  326     _-  36     .-  196    )/   132  5x_2    , 
n-1       \  n-4  n-3  n-2  n-1  n/  1 


and 


6*'   =  (-265     _+  13U6     _   -  385     ,    -706    )/   132 5x2 
n       \  n-3  n-2  n-1  n/  1 
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Relations  for  S1'1  are 


S"i      "  ("  Sl-3  +  6i-2  +  «1-1  -  Si+1  -  Si+2  +  5i+3  )  '   6  SX1   ' 


i  =  3,1*,...  ,  (n-3), 


tit        iii 

6    =0,6    =  0  , 
o         n 


[  5h   -  63  -  252  +  61  +  6q  1 /  6  5xJ   ,  (U.T) 


(55  -6k~   53  +  5o)/  65x?    • 


5  ^  =  6  _  5,  _  5 


1 1  i 
'n-2  ""  I  "  un-5  T  Un-U 


6  ^  =  (  _  5  ^  +  5  ,  +  fi  _  -  6   /   6  6x,   , 

n-j    n  /        1 


and 


i  i  i 


«„-l  ■  "  5„-U  +  S„-3  +  2  6n-2  "  6n-l  "  6n  ) '     6  Sxl 
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Use  of  these  equations  for  finding  5' ,  5"  and  5' ' '  leads  to  mild 

fluctuations  in  <5-values  after  a  few  steps  in  the  Z-direction.   Their  use  had, 

therefore,  to  be  abandoned.   Less  accurate  relations  involving  an  error  of 

0  6^  also  yielded  similar  results. 
1 

v)   The  best  method  for  finding  the  derivatives  of  <5  for  this  problem 
turned  out  to  be  the  use  of  central  difference  formulae 

5!   =  (6i+1-5i_1)/  2  6xx   ,  i  =  l,2,...,(n-l) 

(U.8) 


i    » 

5=6=0 
o    n 


for  the  first  derivative,  and 


5 


i  =  (  6i+l  ~  2  6i  +  5i-l)  '      5xl   •   i=l>2,...,(n-l)   , 


5o  =  (6l-5o   /  2  S4        -  ik'9) 


!! 


6      '  6  _  -  6  )/  2  5x2 
n  =  I  n-1    n  /       1 


for  the  second  derivative.   The  third  derivative  was  found  by  applying  (U.8) 
to  5"  calculated  from  (U.9)«   This  method  gave  the  best  results  for 
5-distributions ,  and  was  therefore  adopted. 
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5.   RESULTS  AND  DISCUSSION 

The  complete  computer  program  was  developed  in  two  stages.   In  the  first 
stage  the  solution  of  eqn.  (3.1)  alone  was  attempted.   This  program  was  tested 
for  the  vapor  and  tube  properties  for  which  results  are  available  in  [21 ] .  As 
pointed  out  in  "Introduction",  Stack  and  Merkle  solve  a  much  simplified  form 
of  (3«l)  since  besides  other  simplifications,  they  use  Cartesian  coordinates 
rather  than  curvilinear  coordinates .   In  fact ,  we  had  modified  our  earlier 
program  to  solve  Stack  and  Merkle fs  equation  (lU)  for  the  condensate  film 
thickness,  and  obtained  identical  results.  The  finite  element  program  for  the 
solution  of  eqn.  (J+.l)  was  then  developed,  and  the  two  were  coupled. 

In  the  following  we  present  results  for  one  vapor-tube  combination  for 
which  experimental  data  is  also  available  [lO,ll]  for  comparison.  Since  our 
analysis  assumes  the  radius  of  curvature  of  the  fluted  tube  to  be  large 
compared  to  the  film  thickness,  it  is  not  realistic  to  consider  sinusoidal 
flutes  with  amplitude-to-pitch  ratio  much  above  0.2.  With  this  in  mind,  we 
computed  results  for  tube  'F'  of  Ref.  [l0,ll].   This  tube  has  seven  skirts  and 
an  I.D.  =  22.9  mm.   Details  of  the  actual  shape  of  the  flute,  and  its 
dimensions  are  not  available  in  [l0,ll],  and  were  therefore  approximated  from 
a  blow-up  of  the  tube  *F'  photograph  in  [ll] .  The  relevant  fluid  (R-113)  and 
tube  (aluminum)  properties  are 


p  =  IU98.3U3  kg/m3 
v  =  3.206Tx10~7  m2/s 
hfg  =  1U5225.56  J/kg 
kf  =  0. 06951  W/m-K 


Twoo  =  318.5  K 


Pv  = 
a  = 

Ts  = 
kv  = 

T„  = 


8.58628  kg/m3  , 

0.01U32  N/m    , 

325.5  K 

205  W/m-K 

318.5  K  (in),  318.8  K  (out) 
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L  =  0.1U2875  m  (  =  -^^  x  |^  )  ,  h  =  0.880TT  mm   , 

P  =  1.6lUl  mm     ,     h0  =  0.155^3  mm    , 

where  kf  and  kw  are  the  thermal  conductivities  of  the  condensate  and  fluted 
tube  respectively,  and  h^  and  h0  are  associated  with  the  flute  shape  (Fig.  3) 
taken  to  "be 


f(x)  =  h.  +  h  cos  ^  .  (5.1) 

to       P 


Thus  hQ  is  half  the  amplitude  of  the  flute.   The  fluid  properties  were 
calculated  at  322  K  (mean  of  Ts  and  Twoo)  from  known  correlations  [25]  .   These 
values  lead  to 


5y  =  0.08051  mm    ,   xp  =  O.876U87 


mm 


After  some  numerical  experimentation  involving  different  step  sizes, 
etc.,  the  solution  of  eqn.  (3.1)  was  started  from  ZQ  =  5xl0~6  with  0  =  1/2 
(corresponding  to  central  differencing  in  X),  5X  =  0.05,  an  initial  <5Z  = 
5xlO~5 ,  and  an  initial  A  =  0.5  (see  eqn.  (3.5)) •   Further  iteration  for  the 
solution  of  (3.1)  was  terminated  when 


^r  >Kk 


<    e  for  all  i  at  every  k  ,        (5«2) 


where  e  was  taken  to  be  10~6.  As  we  marched  downstream  in  the  Z-direction,  6Z 
was  increased  and  the  relaxation  factor  X  had  to  be  decreased  according  to 
Table  1,  which  also  gives  an  idea  of  the  number  of  iterations  before  (5*2)  was 
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satisfied.   Clearly,  it  takes  longer  to  get  the  converged  solution  as  we 
proceed  downstream  in  the  Z-direction.  Reasons  for  this  will  be  apparent 


soon. 


Table  1 Some  Parameters  for  Solution  of  (3.1) 

5Z  X  Z         #  of  iterations 

~  15 
~  30 
~  50 
~  90 


5xlO-6 

0.5 

up  to 

0.002 

8xl0~6 

0.2 

up  to 

0.01 

1.25xl0~5 

0.1 

up  to 

0.0225 

2xl0-5 

0.06 

up  to 

0.086 

Three  divisions  of  the  region  0ABC0  (Fig.  3)  into  finite  elements  for  the 
solution  of  eqn.  (U.l)  were  tried.  The  one  selected  on  the  basis  of  adequate 
computational  accuracy  and  relatively  economical  calculation  had  15  *+ 
triangular  elements  with  108  nodes.  There  were  11  nodes  each  on  the  faces  AB, 
BC,  OC  and  DE  (Fig.  3),  and  15  nodes  on  the  face  A0. 

Fig.  9  shows,  to  true  scale,  the  flute  shape  (thick  curve),  and  the 
condensate  film  shape  (dashed  curves)  at  four  values  of  Z  =  0.001,  0.01,  0.03, 
and  0.086.  As  expected  the  film  thickens  quite  rapidly  in  the  trough  region 
while  it  remains  thin  over  the  crest.  Unfortunately  we  are  unable  to  present 
any  results  for  Z  much  greater  than  0.086  since  convergence  of  our  solution 
for  the  film  thickness  (with  the  same  tolerance  e)  becomes  very  slow.   The 
principal  reason  for  this  slow  convergence  is  the  thickening  of  the  film  in 
the  trough.  Moreover,  it  is  really  futile  to  try  to  speed-up  the  convergence 


is 


since  our  formulation  for  the  condensate  film  thickness  is  on  shaky  ground  as 
the  film  thickens.  The  non-linear  partial  differential  equation  (2.13)  comes 
from  a  boundary-layer  type  analysis,  i.e.,  we  neglect  all  velocity  gradients 
except  those  in  the  direction  normal  to  the  fluted  surface.  While  this  is 
perfectly  reasonable  over  the  initial  portion  of  the  tube  height  when  the 
condensate  film  is  thin,  it  is  increasingly  questionable  as  the  film  thickens 
in  the  trough  region.   It  appears  that  for  the  fluid-tube  combination 
considered  here,  it  is  improper  to  analyze  the  film  in  the  trough  region  in 
the  boundary-layer-sense  for  Z  >  O.O85.   It  is  also  clear  from  Fig.  10,  that 
shows  the  absolute  value  of  d(l/Rj_)/dX]_  at  various  Z  values,  that  this 
derivative  is  negligible  in  the  trough  region  as  compared  to  its  value  in  the 
crest  region.   It  is  therefore  appropriate  to  neglect  the  surface  tension 
effects  in  the  trough  region  once  the  condensate  film  has  thickened.   Fig.  10 
also  plots  the  absolute  value  of  the  derivative  of  flute  curvature  with  X]_, 
and  even  at  the  low  value  of  Z  =  0.002,  it  is  very  different  from  the 
derivative  of  condensate-vapor-interface-curvature  with  X]_. 

While  these  considerations  will  be  undertaken  in  an  extension  of  this 
project,  it  is  encouraging  to  note  from  Fig.  11  that  a  linear  extrapolation, 
on  a  log-log  basis,  of  our  results  to  date  shows  a  fairly  good  comparison  with 
experimental  data.   Fig.  11  shows  the  heat  load  per  half  flute  (in  W)  as  a 
function  of  Z  on  a  log-log  scale.  The  solid  portion  of  the  straight  line  is 
based  on  computed  results  while  the  dashed  portion  is  the  extrapolation.  The 
lone  circled  point  is  based  on  the  experimental  data  [ll]  .  We  also  checked 
that  heat  transfer  rates  across  the  faces  AB  and  0C  (Fig.  3)  matched  within 
1%,     Theoretically  they  should  be  identical  since  faces  OA  and  BC  are 
insulated. 
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6.   SUGGESTIONS  FOR  FURTHER  WORK 

Besides  the  deficiencies  mentioned  above,  our  present  analysis  also 
neglects  the  transverse  curvature  term,  i.e.,  it  assumes  the  film  thickness  to 
be  small  compared  to  the  radius  of  curvature  of  the  flute.  The  analysis  is 
therefore  restricted,  in  its  present  form,  to  flutes  having  a  small  amplitude 
to  pitch  ratio.  This  is,  however,  a  minor  drawback  since  it  seems 
(analytically  at  least)  that  we  can  overcome  this  deficiency  relatively 
easily.   Its  implementation  in  the  computer  code  is  yet  to  be  done. 

Another  assumption,  common  to  all  analyses  to  date,  and  far  more  serious 
in  terms  of  correct  modelling  of  the  practical  applications,  is  the  complete 
neglect  of  vapor  shear  on  the  interface.   Our  code  does  have  the  provision  for 
studying  its  effect  but  since  we  have  not  analyzed  the  vapor  dynamics  yet,  we 
simply  take  the  vapor  shear  to  be  zero  at  present.   One  should  also  consider 
the  coolant  dynamics  for  a  complete  solution  but  all  this  will  undoubtedly  be 
very  demanding. 
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APPENDIX  A 


CURVATURE  OF  CONDENSATE-VAPOR  INTERFACE 
AND  ITS  VARIATION  ALONG  FLUTE  SHAPE 


Equation  (2.13)  for  the  determination  of  condensate  film  thickness 

6(x  ,z)  involves  the  value  of  pressure  gradient  dp/dx  which,  in  turn, 

depends  upon  d(l/R  )/dx  ,  as  is  clear  from  (2.15).   We  present  below  some 

i    1 

relations  for  calculation  of  this  quantity. 

The  basic  terminology  is  present  in  Fig.  3,  wherein,  f(x)  denotes  the 

flute  shape  and  f-j_(x)  the  condensate-vapor  interface.   The  position  vector  to 

-► 
a  point  on  the  flute  is  rw,  and  that  .to  a  point  on  the  condensate-vapor 

-> 
interface  is  r^_ ,  such  that 

♦    ♦    + 

ri  -  rw  =  nw5   >  (A.l) 

where  nw  is  the  unit  normal  to  the  flute  as  shown  in  Fig.  3.  We  also  let  x-^ 
be  the  curvilinear  coordinate  along  the  condensate-vapor  interface. 
Then  by  definition 

5TE  iT*  (A-2) 

1 

Now  ri  =  (x,  fi),  (A. 3) 

l/2 
x      '2 
and  xs    =  /  (l  +  fi  )   dc   , 

b  (A.U) 

where  prime  denotes  differentiation  with  x.   From  equation  (A.U),  we  get 

dxi  ,2  V2 


dx 


=  (1  +  fi  )  (A.5) 
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From  (A. 3)  then 


A 


dx- 


dr.j_   dx 


dx   dXj_ 

-1/2 
=  (1  +  f[2)  (l,fi  ) 


Thus 


dC 


dxi2    (i+fl'2)2 


(-ft  ,  i) 


(A.6) 


From  equation  (A. 2)  then 


'1      I 


(1     +     fi'      2) 


72 


And 


(A.T) 


d      / 1_\    _   d_  /^_\    dx 


dx-j_      Rjl  dx     Rj_      dX]_ 


sgn(fi    ) 


(1+fi'2) 


72 


t       "2    1 


it  t 

fi     -  3- 


f  •    f  • 


(1+fi'2) 


(l+f'2) 


L/2         ' 


(A.8) 


dx 

]  '2    1/2 

since        ^x     =    (l+f      )  in  analogy  with    (A.5K 


(A.9) 


40 


While  equation  (A. 8)  seems  to  be  in  a  convenient  form  for  calculation  of 

dd/R-j  )/dx  ,  it  is  not  practical  since  f,-  is  known  only  at  a  discrete  set  of 
1    1 

t     if       it  t 
unequally-spaced  x-values .   Thus ,  finding  f  j_  ,  f  j_     and  f ^   is  troublesome ,  to 

say  the  least.   Several  attempts  to  use  equation  (A. 8)  had  to  be  abandoned  for 

this  reason  (see  Sec.  U.l  for  details). 

Another  relation  for  the  curvature  of  the  condensate-vapor  interface 

comes  from  use  of  equation  (A.l),  and  the  fact  that  the  unit  normal  vector  nw 

is  given  by 

♦    /   df   dx   x 

nw  =  l~  j — j  1 —  ' 
w      dx^_   dx-j_ 

-l/2 
=  (l+f'2)     (-f\  1)  .  (A. 10) 

Thus,  from  (A.l) 


+    /       6f 


x  -  rfl    »  f  +  1/2  /  '  slnce  rw  =  (x>f > 


(1+f'2)  (1+f'2) 


Therefore,      Jl^/l-fi'f'- ^J-   ,  f'+  6'  .    Sf'f"j/2  j  ,   (A.ll) 

(1+f'2)  (1+f'2)     ' 


i 


where  6  =  85 /3x  but  f  =  df/dx,  etc..   The  reason  for  keeping  3<5/3x  in  place 
of  3S/3x  is  that  solution  of  equation  (2.13)  leads  to  values  of  5  at  equidis- 
tant x  -values  but  at  non-equidistant  x-values.   It  is  therefore  convenient 
1 

to  find  numerically  derivatives  of  6  with  respect  to  x  . 
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Carrying  on  further,  we  can  find  d2rj_/dx2  from  (A. 11),  and  then  use  the 
relation 


1      .    l|ri    I 


^iT  "    (?i'.?i")2    1 


*    ' *   "\P    l  1/2 


R: 


r^ 


(A. 12) 


T  T 

where   r-j_      =  drj_/dx,   to   get   after   some   lengthy  algebra  the   expression 


1  |t56"   +   TU65'    +   y35  +   y252   +   y1    (1-65"   +   26  2    ) 


[yAl+6  2)   +  yv52  +  YQ6] 


Y6 


!/2 


7  T8 


(A. 13) 


where  y     =  f 


Y2  =    (— V) 

1+f  2 


Y3   = 


-2f 


"2 


(1+f   2) 


72 


-V2  1     "o 


1+f 


Y5    =    (1+f   2) 


Yg   =    1+f   2       , 


•"    \2 


1+f 


2f 

y8  = —r 


(l+f  2) 


/2 
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It  is  interesting  to  note  that  many  theoretical  analyses  that  clain  to 
determine  the  surface  tension  effect  'based  upon  the  actual  shape  of  the 
condensate-vapor  interface  invariably  use 

t» 
1       5 


*1  =  (1+6'V/2 


(A.lU) 


for  the  curvature  of  the  interface.   This  follows  from  the  correct  relation 

(A. 13)  if  only  the  first  terms  in  both  the  numerator  and  the  denominator  are 

retained.   However,  even  for  small  6,  there  is  hardly  any  justification  for 

dropping  all  other  terms  in  (A.13).   The  difference  between  the  actual 

d(l/R-,-)/dx  and  that  obtained  from  (A.lU)  gets  further  amplified  due  to 
1 

differentiation . 

For  convenience  let  us  write 


1  =  jNl    N.sgn(N) 


where  the  expressions  for  N  and  D  are  clear  from  equation  (A.13).   It  is  then 
easy  to  get 


d  ,1  v  _  sgn(N) 
dx  ^V      *72~ 
1        D 


6"'(y5-6y1)  +  6"{y5'  +  5(yu-y1')  +  3y^  }   +  Y1' 


+    6   2(yu+   2Yl')    +    6'{Y3+    5(2y2+    Yu  '  )  }  +    <52Y2'    +    5y3' 


2^-  {6'(2y65"  +   2y?5  +   yQ)      +      (l  +    S^hg'    +    52YT'    +    5y8'},    (A.15) 


A3 


Where 


dYi  ,    -1/2 

Y,       =  1 —     =   f        l+f   2 
1  dx]_ 


dY2  3f"2 


'2  dx. 


ni  2f   f  2 


1        (1+fV/2L  1+f2j 


dY3 


2f 


3  dx.  ,      2 

1  (1+f  2) 


3f  f 
1+f  <■ 


-  2f 


dY, 


•IV 


h     "  dx 


i„    2 


1  1+f  2  (l+f'2) 


3f  2   +     Yf  f      )    +     -2 - 


2f"3 


(1+f'2)3 


dYc 
Y5     =  dxf  "  3f  f 


dY 


r6     "  dx. 


6   =  2f'f"(l+f'2) 


-V2 

TO  '    *■ 


dY 


T  2f 


'T  dx. 


(1+f  ') 


it  i    itp 

,    »t      2f  f        v 

5/2    lf      "  ^ ' 

'2*    '  1+f 


and 


r8   _  dZ 


f  f 


'2    v  '2 

1+f  1+f 


-   f      ) 
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Iquation  (A. 15)  is  quite  complicated  but  as  mentioned  in  the  discussion  on 

3  i   3 
qn.  (2.13),  it  contains  3  6/3x]_  .   It  is  therefore  essential  to  determine  the 

ondensate  film  thickness  as  accurately  as  possible.   Otherwise,  errors  in 

inding  the  third  derivative  of  5  w.y«t.  X]_  numerically  will  be  excessive. 

ue  to  these  numerical  difficulties,  it  was  found  that  use  of  eqn.  (A. 15)  to 

ind  d(l/Rj_)/dX]_  resulted  in  mild  oscillations  in  5-values  after  a  few  steps 

n  the  Z-direction.   Several  methods  (at  least  five)  were  used  to  find 

-derivatives  w.r.t.  X]_  but  all  resulted  in  oscillatory  6-values;  the  least 

scillations  resulted  from  the  use  of  central-difference  formulae. 

In  order  to  overcome  this  difficulty,  another  relation  for  1/Rj_  and  its 

erivative  w.r.t.  X]_  was  developed.   This  relation,  based  on  eqns.  (A. 2)  and 

A.l),  is 

32(rw  +  n^S) 


Ri 


3x, 


32rT,    tfn 


■  v 


+  fi- 


v 


3n 


»*i2 


»«i2 


+  2- 


v 


3x. 


36    -v   326 


(A. 16) 


here  3 /3x-^  has  been  approximated  by  3/3x]_.   Noting  that  rw  =  (x,f),  and  n^. 
3  given  by  eqn.  (A.10),  we  get  from  (A.l6),  after  some  algebra,  that 


(ax+  a25  +  a3S'  +  0^5  r  +  ( 01  +  ?>2  5  +  0  6'  +  0^5") 


»\2 


7  2 


(A. IT) 
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where   <5      =    85/3x        ,      f     =  df/dx,    etc.,    and 
1 


a     =     -3  f'      , 

1  1 


a2    =  - 

a3   =     -23. 


1 Aff2  *"M 

'ox/2        1   +  f  2 


i 


aU   =     -B.    f        , 


0!    =   f    (1+f   2) 


B2   =  o2f     -   BiBl;f 


83   =  2ax      , 


and        .  , ,       'of    /2 

3U  =  (l+f  d) 


From  equation    (A. IT)    it    follows   that 
d    /    1 


K  -(p— )    =  Ri{(o1  +   a26  +   a36'    +   au6")(S2  +   c^'s  +  Sc^a'    -   3816"  +    c^fi"') 
+  (B1  +  B26  +  B36'   +   Bu6")(-a2  +   B2'6  +  382<s'   +  3(^5"  +    8^  s"' )  }      ,      (A.l8) 


J  A 


where 

dot  "  /  i      "o  t     It  ? .                               n     '  o    "o 

«    =  uct2  =       1                f  (Uf  2  +  13f  f      )        fiv  _  28 f  2f  3 

^       '   dXl  ,          'ov3                    1  +   f'2                                   ,         '0    2 


(1+f  2)  (1+f  2)' 


and 

dp  2 

'  O  till,  .  ^       T»  » 

3       =       2     -  a     f    +  f  (23,  a    -  a  3  )  -  6  3,    f 

2  ^ —  2  l»   2  1   1  1   h 

1 

3  3 

Equation    (A.l8)    also  contains    3    6/8x      ,   "but  unlike   equation    (A. 15),    it   yields 

1 

stable  (non-oscillatory)  values  of  6  upon  integration  of  eqn.(2.13).   This 
relation  was  therefore  used  in  the  final  computer  code. 
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APPENDIX  B 


Nonlinear  Equations  Solution  by  Newton-Raphson  Method 


Consider  the  set  of  equations 

[K]{x}   =    {R} 


n 

k^  ( x 

1'      2'  ll     J  x     1"      2 


or        Z      kij(x    >    x    ^    ••>    xn^xj    =   Ri^xn  '    x^ »    ••»    xn^      '    i  =   -1'   2'    "»   n" 


j=l 


where  k,-  <    and  R,-    are   non-linear   functions   of   x    ,    x    ,    ..,    x_ . 
j-j  j-  IP 

n 

Let   F,-(x    ,   x  <,..,   xn)   =     E     k-n(x    >   x    >    •  •  »   Oxi    -  Ri^x    >   x    >    ••»   xn  K 
i12  n  .    1     X«J      1        2  nj  x      1        2  n 

i   =   1 ,   2,    . . ,    n. 


th 
So  Fj_   is  the  residue  in  the  i   eqn.   For  solution  Fj_  =  0,  i  =  1,  2,  ..,  n. 

Expanding  Fj_  by  Taylor  series ,  we  write 

F,-(x  +  Ax  ,  x  +  Ax  ,  . .  ,  xn  +  Axn ) 
1   1     12     2       n     n 


=  Fj_(x  ,  x  ,  . .,  xn)  +  E   dFi  (x  ,  x  ,  ..,x   )Ax<  +  ..., 


12       "    <=i  ^77   1   2 


j=l  3xj    ■*■  n 

i  =  1,  2,  ..,  n. 

Neglecting  all  higher-order  terms  and  setting  leftside  =  0  for  a  solution, 

we  get 

n 

£   <*Fi  (x  ,  x  ,  ..,  xn)Axi  =  -F,-(x  ,  x  ...,  x.)  ,   i  =  1,  2,  ..,  n. 
j=13xT   1   2  J        1   2 
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Thus  Newton-Raphson  iteration  algorithm  can  be  written  as 

,  ,m    m+1         m 
fj]  {Ax}     =  -{F} 

m+1       m       m+1 
{x}      =  {x}   +  {Ax} 

where  m  denotes  m^"  iteration  &  [j]   is  the  system  Jacobian  at  the  m^"1  itera- 
tion. 

r  im 
The  elements  of  I J]   are 

3F- 

Tm   _    1  /  m     m        m\ 

ij  ~      dx-         1    '  *2   *  ""  '   "   " 

m 
If  it  is  too  time  consuming  to  compute  [J]   at  every  iteration,  find  [j]  once 

and  use  it  until  convergence,  which  will  be  slow,  however. 
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TUBE  SUPPORT/CONDENSATE/4!  JjJ 
DRAIN-OFF  PLATES 

ii  i*; 

AIR  COLLECTION  DUCT -U' 

STEAM 
INLET 


VnrrrTii^p 


COOLING 
WATER  INLET 

HOTWELL 


RETURN  HEADER 


VERTICAL  FLUTED  TUBES 


AIR  REMOVAL 

INE 

DOWNCOMERS 


COOLING  WATER 
OUTLET 


CONDENSATE 
OUTLET 
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Fig.  1.  Vertical  Submarine  Condenser 
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Fig.  4    Control  Volume  in  the  Condensate  Film. 
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Fig.  II.    Heat  Load(inW)    as  a  Function  of  Z  .—, Computed  j 
Extrapolated-,©, Experimental  [ll_ 
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